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ABSTRACT 

Two charge- and current neutral plasma beams are modelled with a one-dimensional PIC 
simulation. The beams are uniform and unbounded. The relative speed between both beams 
is 0.4c. One beam is composed of electrons and protons and one out of protons and negatively 
charged oxygen (dust). All species have the temperature 9 keV. A Buneman instability develops 
between the electrons of the first beam and the protons of the second beam. The wave traps 
the electrons, which form plasmons. The plasmons couple energy into the ion acoustic waves, 
which trap the protons of the second beam. A proton phase space hole grows, which develops 
through its interaction with the oxygen and the heated electrons into a rarefaction pulse. This 
pulse drives a strong ion acoustic double layer, which accelerates a beam of electrons to about 
50 MeV, which is comparable to the proton kinetic energy. The proton distribution eventually 
evolves into an electrostatic shock. Beams of charged particles moving at such speeds may occur 
in the foreshock of supernova remnant shocks. This double layer is thus potentially relevant for 
the electron acceleration (injection) into the diffusive shock acceleration by supernova remnants 
shocks. 



Subject headings: physical data and processes: acceleration of particles 
plasmas — methods: numerical 



physical data and processes: 



1. Introduction 

The blast shell of a supernova remnant (SNR) 
can accelerate charged particles to cosmic ray 
energies through diffusive shock acceleration by 
the SNR shocks (lUchivama et al1l2007t iRevnoldsl 
2008; Hillasl 2005 ). Particles are scattered up- 
stream and downstream of a shock and some 
particles cross the sho c k repeatedly, gaining 
energy each time (iBeil 1 19781; iKrvmskvl Il977t 
Blandford fc Ostrikerll97alMalkov fc Drurvh oOl) 



Diffusive shock acceleration requires a seed pop- 
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ulation of particles with energies well in excess 
of the thermal ones, since only such particles can 
cross the shock repeatedly and be scattered ef- 
ficiently by the magnetohydrodynamic waves on 
either side of the shock. The presence of shock- 
accelerated electrons is evidenced by the radio 
emissions of SNRs. It is not clear yet how the 
electrons can reach the energy threshold for dif- 
fusive shock acceleration; that requires that their 
kinetic energies are comparable to those of the 
ions. Some electrons must be accelerated locally, 
that is close t o the SNR shock or close to one of 
its precursors (|Revnolddl200i lHilladl2005l ). 

This injection of electrons is thought to be ac- 
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complished by the i nterplay of energeti c ions with 
electrons and dust ( Ellison et al. 1997h bv means 
of pl asma instabilities (jCargill fc Papadopouloj 
1988). The instabilities occur in the upstream 
plasma, just ahead of the SNR shock. This 
region is called the foreshock. One source of 
the ion beams is the downstream plasma. Ions 
can leak through the shock boundary and out- 
run it (|Malkovl Il998h . The shock can also 
reflect a substantial fraction of the upstream 
ions. If the reflection is specular, then the 
ions can cross the foreshock at twice the shock 
speed. Reflected and leaked ion beams can 
be observed in situ at the Earth bow shock 
(|Eastwood et al.ll2005l ). Particlc-in-cell (PIC) sim- 
ulations of fast plasma shocks suggest, that they 
may also occur at SNR shocks when the upstream 
magn etic field is orthogo n al to the shock nor- 
mal ( Hoshino k. Shimadal 2002]: Lembege et al 



2004 t IChapman et all 12003" lAmano fc Hoshino 
2007). Precursor shocks move ahead of the main 
SNR shock. They can reach mildly relativistic 
speeds for the most e xtreme supernova explosions 
(Kulkarni et al Such shocks can also reflect 

ions i n the presence of a quasi-parallel magnetic 
field (jDieckmann. Shukla fc Drury||2008h . 

We focus in this work on ion beam insta- 
bilities in the foreshock, but we do not model 
the shock itself. We rely on a foreshock model 
that has been widely u s ed, e. g. in the works 
bv IShimada fc Hoshinol d2003l . |2004 and by 
Dieckmann et al.l (|2007l l2008al) . This model con- 
siders an ion beam, which moves through the up- 
stream plasma, consisting of cool electrons and 
protons. The system is spatially uniform. Of- 
ten a second, counter-propagating ion beam is 
introduced, which cancels the current of the first 
one. This second beam can be due to ions that 
have been reflected by a perpendicular shock 
and which return to the shock, after they have 
been rotated by the u pstream magnetic field 
( M c Clements et al 1 1l997h . Here a beam of co- 
moving negatively charged oxygen ions cancels 
the current of a proton beam. 

This foreshock model can be used if all pro- 
cesses develop on scales well below an ion Lar- 
mor radius and an inverse ion cyclotron frequency 
in the case of a perpendicular shock. The fore- 
shock model may be representative over larger 
scales otherwise, because ions can expand much 



farther if they move along the magnetic field, as 
it is observed at the Earth's bow shock. This 
foreshock model is not self-consistent, because it 
decouples the foreshock physics from that of the 
shock. It also assumes a plasma that is spatially 
uniform, while the shock reflected ion beams vary 
in their mean speed and temperature. However, 
we can examine with this model the acceleration 
of the upstream electrons as a function of the bulk 
parameters, such as the beam temperatures and 
mean speeds and the background magnetic field. 
We can also separate individual processes from 
the global shock d ynamics, that involves coupling 
over many scales ( Hoshino! 120031 ). This simplifies 
their identification and interpretation. The peri- 
odic boundary conditions also allow us to reduce 
the box size, making PIC simulations computa- 
tionally efficient even for realistic ion-to-electron 
mass ratios. The insight gained from such simula- 
tions can then guide us in the selection of plasma 
parameters for the more realistic and demanding 
PIC simulations of full shocks. 

Simulation studies employing this foreshock 
model have revealed that proton beam-driven in- 
stabilites are not efficient electron accelerators. 
An obstacle is the mass difference between the 
electrons and the protons. Most plasma instabili- 
ties saturate by accelerating the electrons. The 
fields are too weak to accelerate the ions and 
the heated electrons quench most of the known 
plasma instabilities, e.g. through Landau damp- 
ing. Often only a few per cen t of the ion energy is 
trans f ered t o the electrons (IShimada fc Hoshinol 
120031 120041 : iDieckmann et all 120071 l2008al ). Ini- 
tially, this is also the case for i on beams that ex 



pand i nto the upstream plasma (|Dieckmann et al 



2008b) prior to the formation of a shock. Certain 
mechanisms dKuramitsu fc Krasnoselskikb 2005: 
Katsouleas fc Dawsonlll983l ) can, in principle, ac- 
celerate the electrons to relativistic speeds, but 
this remains yet to be demonstrated with self- 
consistent PIC simulations. However, the elec- 
trons, that have been accelerated by a n instabil- 
ity li ke the Buneman instability (BI) (jBuneman 



19581 ) , are often spatially b unched, for example in 



form of phase space ho le s dRoberts fc Berkl 1967 ; 
Korn fc Schamellll996al rbl: IDieckmann et al.ll2007[) . 



Electron bunches or plasmons have a much higher 
inertia than electrons and they modulate the ions 
to give ion acoustic waves. The ion acoustic waves 
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can trigger further processes, one of which we 
consider in this work. 

We investigate a plasma consisting of two spa- 
tially uniform beams. One beam is composed 
of protons and electrons and the second one of 
protons and oxygen carrying a single negative 
charge. Each of the beams is charge and cur- 
rent neutral. Both beams move at a relative 
speed of 0.4c. All species have the temperature 
9 keV. These initial conditions are the simplest 
ones that take into account three particle species 
and they permit a first assessment of the impact 
of heavier ions. Section 2 discusses in more detail 
these initial conditions and their potential rele- 
vance for a foreshock plasma. Section 3 describes 
the simulation results, which can be summarized 
as follows. A BI develops between the electrons 
and the protons of the second beam. The BI 
saturates by forming electron phase space holes 
and these plasmon s feed wave energy into the 
ion acoustic wave (|Mendonca fc Bingham 2002; 
Mendonca et al.l 120051 ). Proton phase space holes 
develop, once the ion acoustic waves are strong 
enough to trap the beam protons. The electrons 
eventually thermalize and they reach a mildly rela- 
tivistic thermal speed, which is stabilizing the sys- 
tem for some time. The BI plays no further role 
after the electron thermalization. We examine in 
detail a proton phase space hole that grows into 
bein g a proton density pulse with a negative den - 
sity llnfeld et alJll989t llnfeld fc Rowlands! !l990). 



The growth of this proton phase space hole may 
be c aused by its interaction w ith the negative oxy- 
gen (jEliasson fe Shuklall2004l) or by th e disruption 
of the electron current ( Smith! 19821 ). The pulse 
triggers a strong ion acoustic double layer and it 
evolves subsequently into an electrostatic shock. 
The electron phase space structure and evolution 
of this double layer rese mbles that, found by a 
recent simulation study bv I Newman et al. ( 2001 ). 
That study addressed the electron acceleration in 
the Earth's auroral ionosphere. However, our sim- 
ulation yields the much stronger electron accel- 
eration that is necessary for injecting the elec- 
trons into the diffusive shock acceleration at SNR 
shocks. The proton phase space distribution also 
differs from the standard picture of double lay- 
ers outlined, for exa mple, in the review papers b y 
Smith! (|l982h and bv lRaadu fc Rasmussenl (|l988l ). 
The trapped and free streaming protons in our 



simulation are well-separated in velocity, which is 
not generally the case for double layers. The sub- 
stantial free energy stored in the difference of the 
mean velocities of both beams causes the extreme 
electron acceleration up to 50 MeV. We discuss 
our simulation results in section 4, we relate it 
to previous investigations of double layers and we 
discuss the necessary future research. 

2. The initial conditions and the simula- 
tion method 

The bulk plasma of a supernova blast shell can 
expand into the ambient medium at a speed of 
up to 0.1c-0.2c, at least during its ea rly phase 



(|Vink et al.l 120061 : iFransson et all [2002! ). A frac 



tion of the blast shell plasma may expand even 
faster. A plasma flow speed of 0.9c has, for 
example, been observed for a subshell ejected 
by a particularly y iolcnt supernova explosion 
(Kul karni et al. I ll998h . Shocks form where the ex- 
panding plasma impacts on the ambient plasma, 
the interstellar medium. The typical shock speeds 
should be comparable to the few per cent of c the 
main blast shell plasma has, but some precursor 
shocks may be considerably faster. The detailed 
structure in particular of the fastest SNR shocks 
is unclear. We assume here, that the foreshock 
region is resembling at least qual itatively that of 



the well-known Earth bow shock ( Eastwood et al 



20051 ) and we focus on the ion foreshock. The ion 
foreshock is characterized by the presence of ion 
beams that can propagate far into the upstream. 

However, we expect that a hot background of 
electrons is present in the foreshock, as PIC sim- 
ulations of mildly relativistic shock s demonstrate 
(|Dieckmann. Shukla fc Drurvll2008l ). 

We examine the interaction between two spa- 
tially uniform and unbounded plasma beams, 
consisting of one negatively and one positively 
charged species each. One beam consists of elec- 
trons and protons and the second one of protons 
and oxygen carrying a single negative charge. 

The simulation box frame is the reference 
frame, in which both plasma beams move in op- 
posite x-dircctions and with equal speed moduli 
Vb- The simulation models only the x-direction, 
which allows us to resolve simultaneously the ion 
scales and the electron scales at the correct mass 
ratio. 
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The beam 1 consists of electrons with the 
plasma frequency Q e — (e 2 n e /m e eo) 1 ^ 2 , where e 
and n e are the elementary charge and the electron 
number density and eo is the dielectric constant. 
The equally dense protons have the plasma fre- 
quency il p = Rp 1 ^ 2 £l e , where R p — m e /m p . Both 
species move at the mean speed < v >= —Vb 
with vi, — c/5. The beam 2 is composed of pro- 
tons with £l p and oxygen, which carries a single 
negative charge. The oxygen plasma frequency is 
fl = Qp/4. The mean speed of both species of 
beam 2 is < v >— Vb- All species have the tem- 
perature T = 9.1 keV, which gives the electron 
thermal speed v e = (fcsT/m e ) ' 5 = Vb/1.5. By 
neglecting the T ^ and the heavy species of 
each beam, the dispersion relation for the beam- 
aligned electrostatic modes can be cast into the 
form 



1 



<> 2 



n 2 



(uj + kvb) (u — kvbY 



= 0. 



(1) 



The wavenumber and frequency in the reference 
frame of beam 1 of the most unstable wave are 
k u = Q e /2vb and lo u = £! e . The Doppler-shift will 
reduce the frequency of the most unstable wave to 
f2 e /2. The wavelength A u = 2%/k u . The phase 
speed v p h = oj u /k u w 2vb in the reference frame of 
beam 1 or Vb in the box frame. Replacing the oxy- 
gen with electrons gives symmetric beams, which 
lets the phase speed of the waves vanish in the box 
frame. An exact resonance between the plasmons 
and the proton beam would not be possible. The 
estimates for ui u and k u , as well as the growth rate 

n Rn /n F ~ (3 1 / 2 /2y 3 )(n n /n„) 2/3 « 0.056 of the 

cold BI (jBuneman 19581 ) are fairly accurate. Fig- 
ure [1] shows the numerically obtained growth rate 
£Ib /f2 e f° r our plasma parameters. This growth 
rate m ap has been obtai ned by using a warm fluid 
model (jBret et al.l 120081 ). The growth rate map is 
symmetric around w&, so we do not need to con- 
sider k y . The growth rate in the warm plasma 
model is lower than in the cold one by about 10%. 
The corresponding wave number is increased by 
20%. The growth rate decreases for increasing k z . 

The PIC simulation method Dawson (|l983l) 
solves the Vlasov-Maxwell set of equations with 
the method of characteristics. It approximates 
the phase space density /(x, v, t) by an ensemble 
of computational particles (CPs). The CPs corre- 
spond to phase space volume elements and are not 



physical particles. The charge q cp and mass m cp 
of a CP can thus differ from those of the particles 
it represents. The charge-to-mass ratio q C p/m C p 
must, however, equal that of the physical particles 
and the ensemble properties of the CPs approxi- 
mate well those of the physical plasma. The tra- 
jectories of the CPs are prescribed by the Lorentz 
equation of motion. The Maxwell-Lorentz set of 
equations can be normalized, which provides us 
with results that are independent of the plasma 
density. 

Let e and m e be the elementary charge and the 
electron mass. 

Our code is based o n the numerical scheme 
outlined by lEastwood (1991) and it solves the 
Maxwell's equations for all components of the elec- 
tromagnetic fields together with the Lorentz equa- 
tion, which updates the three-dimensional vector 
of the relativistic momentum for an ensemble of 
computational particles (CPs). 

The relevance of double layers for solar and 

astrophysical particle acceleration has been pro- 

posed first bv lAlfven fc Carlauistl (|l967l) . iLembege fc Dawsonl 
(1989) have found double layers with PIC simu- 
lations of oblique shocks, which dissipated the 
electron energy. In our PIC simulation the double 
layer emerges out of the BI and prior to the for- 
mation of the shock. A relativistic and monoener- 
getic electron beam is ejected by the double layer. 
The time evolution of the proton pulse implies, 
that the mean energy of this beam rises in time to 
about 50 MeV, which is a much stronger accelera- 
tion than that reported, for examp le, for previous 
PIC simulations of double layers (jSato fc Okuda 
19801) . 



The nonrelativistic beam speed and the high 
temperature ensure that the BI grows fastest. The 
approximation of this system by a ID PIC simu- 
lation, which enforces a monodirectional (plane) 
wave spectrum, may thus initially be justified. 
This geometry and the idealized initial conditions 
facilitate the interpretation of the physical pro- 
cesses and they permit us to resolve the relevant 
spatio-temporal scales. However, the growth rate 
of the oblique modes is not negligible and the wave 
power is isotropized in tim e through weak turbu- 



lence ( Ziebell et all 12008 ). The PIC simulation 



may also become unrealistic when the electrons 
have been accelerated to relativistic speeds. This 
is, because of the unresolved higher-dimensional 
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instabilities, the generation of energetic photons 
through synchrotron- and bremsstrahlung and the 
back-reaction of these processes onto the plasma 
( Weibellll959l:lBret A.ll20o5 ISchlickeiser fe Lerche 



2007 ~ ISchlickeiserl 120031; iFleishman fc Toptyg 
20071) . 

Charged dust and heavy ions exist in the su- 
pernova ejecta and in the interstel lar medium 
(jEllison et aIlll997l : lMeikle et al.ll2006l) . Their ion- 
ization state and their relative densities are un- 
clear. It is typically only the dust, which has 



a negative time-dependent ch arge (jMomeni et al 
20071 IVerheest fc Pillay[ l2008h . while the oxygen 



20081 ). Our 



is positively charged (|Katsuda et al. 
initial conditions are thus not necessarily rep- 
resentative for SNR blast shell plasmas. How- 
ever, the purpose o f the oxygen is to intro duce 



a slower time scale (jEliasson fc Shu kla 2004) and 
also to obtain a resonance between plasmons and 
the io n acoustic wave, which facilit ates their cou- 
pling (jMendonca fc Binghamll2002r ). Both condi- 
tions can probably also be achieved by heavier 
ions, by a mix of ions or by dust. 

The beam speed of 0.4c is too high for typi- 
cal SNR bulk flows, which ca n reach about 10- 
20% of the speed of light c (jVink et all 120061: 
Fransson et al.l2002l ). However, the shock-reflected 



ion beams can outrun the shock. Some compo- 
nents of the expanding plasma may also move 
faster than the mean speed of the SNR ejecta. 
Precursor shocks of particul arly strong SNRs ma y 



reach speeds of up to 0.9c (jKulkarni et al.lll998l 



We select this beam speed, because it facilitates 
the processes we consider and because it allows 
us to estimate the maximum energy, which the 
electrons could gain in typical SNR flows through 
the plasma instabilities discussed here. The beam 
speed is also low enough to favor the growth of 
the BI rather than that of the multi-dimens ional 
filamentation instability ( Lee fc Lampei 19731 ). 

The simulation box length is L m 1170A U . It is 
resolved by N g = 7 ■ 10 4 grid cells, each with the 
length ss 0.31Ad. The electron Debye length 
is Ad = v e /Q e . Each of the four species is repre- 
sented by 250 CPs per grid cell, giving a total of 
N p = 1.75 • 10 7 per species. The system is followed 
in time for tsile w 2600, which is subdivided into 
10 5 time steps A f . The beams propagate the dis- 
tance tsvi, < L/5 and even a light pulse can cross 
L just once. Effects due to the periodic boundary 



conditions of the simulation box can not occur. 
3. Simulation results 

The electrostatic BI develops between the 
electrons of beam 1 and the protons of beam 
2. However, the secondary instabilities give 
rise to electric fields, which will let all par- 
ticle species interact. Figure [2] displays the 
time evolution of the relevant energy densi- 
ties, which are the energy density Ep(t) = 
eoNg l Y!jl\El{jA x ,t) of the electrostatic field 

and the E t (t) = iV^A" 3 rm^^Vj) - 1], 
where the summation is over all particles j of the 
species i with the mass mi and a speed Vj (t) giving 

T( Vj ) = (l-v]/c 2 )- 1/2 . 

The Ep oc E% grows in the interval tfl e < 100 
at an exponential rate 2Qi « 0.06, which is below 
2Hb = 0.098. The growth rate l^s is calculated 
for a sine wave. The L » )„ and the broad k x 
spectrum of unstable waves in fig. [T] implies, how- 
ever, that we integrate over a mixture of waves, 
all of which have a growth rate < Hb and the 
initial growth of Ep cannot be compared directly 
to Qb- The electron trapping lets Ep(t) saturate 
at til e w 100 and it causes the growth of E\(t). 
Thereafter Ep(t) decreases and reaches a meta- 
equilibrium at til e w 10 3 . The E\(t) is almost 
constant until t£l e w 1500, when it and Ep(t) start 
to grow again. All ion beams have lost less than 
1% of their energy prior to tVL e w 1500, confirming 
that the BI is not an efficient electron accelerator. 
Then E 2 (t) and E 3 (t) start to change and even 
E±{t) has decreased by 1% at tfl e w 2600. 

Figure [3] shows the phase space distributions of 
the trapped electrons and of the protons of beam 
2 at ttt e = 145 in a subinterval of the simula- 
tion box. The trapped electrons gyrate around a 
potential that is moving with the mean speed u& 
of beam 2, which is also the phase speed of the 
Buneman wave. Only a fraction of the electrons 
is trapped. The electrons of the untrapped bulk 
population move on oscillatory phase space paths. 
The electron phase space hole in the interval < 
x/X u < 1 has the expected size, while the electron 
phase space hole in the interval 2<x/A M <3.5is 
larger. Electron phase space holes increase their 
size by coalescence and their mutual inte raction 
will deform them (jRoberts fc Berkj Il967[ ). The 
mean momentum of the proton beam is weakly 
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modulated. The initially periodic train of elec- 
tron phase space holes will coalesce, until only a 
dissipative equilibri um and solitary electron phase 
space holes remain dKorn fc Schamellll996allbT l. 

Figure [4] confirms this. A simulation box in- 
terval is selected at the time tfl e = 10 3 , in which 
an electron depletion at ;r/A„ « 9 separates the 
turbulent region x/ \ u < 8 from a smooth plateau 
at larger x. The latter constitutes a dissipative 
equilibrium. The supplementary movie 1 follows 
the time-evolution of the electron phase space dis- 
tribution at the position of the electron depletion. 
The color scale is the 10-logarithm of the number 
of CPs. Note that the absolute position in the 
box is given in movie 1, while fig. 2] uses a rela- 
tive position. The movie 1 demonstrates that it 
is this density depletion, which is responsible for 
the turbulent phase space region. A spatial corre- 
lation is visible in the fig. |4] between the electron 
depletion and a phase space hole in the protons 
of beam 2. This explains why the electron deple- 
tion is stationary in the reference frame of movie 
1, which moves with Vb in the box frame. The 
electron temperature is now relativistically high, 
by which the linear damping of the ion acoustic 
waves is reduced and these modes can accumulate 
wave energy. The thermal anisotropy of the elec- 
trons causes the growth of the Weibel instability 
(IWeibellll959h . )ut it is suppressed here by the ID 
geometry. 

Figure [5] plots the E x , the density of the pro- 
ton beam 2 and that of the electrons at the same 
location and time as fig. |U The electric field at 
x/X u « 9 is attractive for the protons and repels 
the electrons, which is in agreement with the ob- 
served phase space structure in fig. [U The electric 
field amplitude at this position is not unusually 
large, while the plasma density depletion is. The 
oxygen density is unaffected at this time, as we 
show below. 

The component energies in fig. [5] demonstrate 
that the equilibrium between Ep(t) and E\(t) 
breaks at t£l e « 1500 and that a secondary in- 
stability sets in after that time. An examination 
of the full phase space data shows, that this sec- 
ondary instability is spatially localized. Two simi- 
lar, but spatially well-separated structures develop 
in the phase space distribution of the protons of 
beam 2 in the box. That two structures rather 
than one grow suggests, that they develop with a 



reasonable probability and that their growth is not 
accidental. Figure [6] displays the larger of these 
two structures at ttt e = 2000. A different sub- 
interval of the box than in the figd]is investigated, 
since the proton phase space structures convect 
with V},. We refer to this proton structure as a 
pulse. This pulse constitutes a density dep letion 



pulse , inis pulse constitu tes a density dep letion 
(Infol d et al. . 1989; llnfeld fc Rowlandsl l l990l ). The 
supplementary movie 2 animates in time the devel- 
opment of the pulse out of the proton phase space 
hole. The color scale shows the 10-logarithm of the 
number of CPs. The movie 2 displays absolute box 
coordinates and its reference frame moves with Vb 
through the box. It demonstrates that the pro- 
ton phase space hole is stationary for a long time, 
before it develops into a rapidly growing pulse. 
It is this pulse that is responsible for the energy 
changes in the fig [2 

The electron phase space distribution reveals 
in fig. [5] a double layer that is driven by the pulse. 
This double layer accelerates on 5X U a beam of 
electrons out of the bulk to about 7.5 MeV. The 
movie 1 demonstrates that the mean energy of the 
electron beam grows in time, while the beam tem- 
perature remains unchanged. The electron beam 
then interacts with the bulk electrons through 
a two-stream instability, which is responsible for 
the large beam momentum oscillations. The elec- 
tric field of the double layer is strong enough to 
decrease the mean speed modulus of the proton 
beam 1, which supplies the energy for the electron 
acceleration. The movie 2 is stopped, when the 
protons of beam 1 and 2 start to mix by the for- 
mation of an electrostatic shock. The movie 1 is 
stopped, when the electron beam energy exceeds 
the one resolved by the movie window. 

The electron distribution at tVl e = 2600 is illus- 
trated by fig [7] The double layer has not spatially 
expanded but the electrons are now accelerated to 
50 MeV in the box frame, exceeding even the 20 
MeV energy of the protons in the same reference 
frame. The total density of the ultrarelativistic 
beam in fig. [TJa) is about 40% of the density of the 
bulk electrons in fig. Efc). The momentum distri- 
bution of the bulk electrons is practically identi- 
cal to that in the figlH This ultrarelativistic elec- 
tron beam would drive the unresol ved and faster - 
growing oblique mode instability (jBret A. 2006), 



rather than the two-stream instability evidenced 
by the movie 1. The type of instability should, 
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however, not be important for the evolution of the 
electron double layer because it forms well behind 
it. 

Figure [5] displays the phase space distributions 
of both proton beams in the same spatial inter- 
val as fig. [7] The proton distributions of both 
beams have merged and a large electrostatic shock 
forms at x/X u » 55. A second electrostatic shock 
involving only the protons of beam 2 occurs at 
x/X u « 92. The movie 2 illustrates that the sec- 
ond shock is driven by the rapid expansion of the 
proton density pulse. It is remarkable that some 
protons of beam 2 have increased their p x momen- 
tum by a factor 2.6 through the electrostatic field, 
which is also shown by fig|5] The electron double 
layer is thus also a proton accelerator. The elec- 
trostatic field of the double layer does not change 
its sign, which contrasts the bipolar fields of phase 
space holes, like in fig [5] 

Solitary phase space holes ar e typically time- 
stationary in one dimensio n (jRoberts fc Berk 
19671 : iKorn fc Schame]|ll996al |bl) and the question 
arises, why the proton phase space holes grow here 
into such powerful rarefaction pulses. The proton 
phase space hole is stable over a long time and it 
then suddenly turns into a pulse, which suggests 
an effect that develops on a slower time scale. The 
electric field of the proton phase space hole will 
gradually displace the oxygen ions. Figure [9] dis- 
plays the density distribution of the oxygen ions 
at the three times tft e — 1000, ifi e = 1600 and 
til e = 2000. The frame is moved with the beam 
2 and the proton phase space hole is located at 
(x — Vbt) w 20 in the chosen frame. The oxygen 
depletion at tfl e — 10 3 is comparable to the sta- 
tistical fluctuations of the number of CPs. The 
integration of the density over 3 grid cells reduces 
these statistical fluctuations by the factor V^j but 
it keeps the density depletion value unchanged. 
The proton phase space hole grows into the pulse 
at tilp sa 1500. Figured) reveals that the oxygen 
density modulation is significant at around this 
time. It is thus likely that the sudden growth of 
the proton phase space hole is a consequence of 
this oxygen depletion. As the proton pulse grows, 
the structure and amplitude of the oxygen de- 
pletion remain qualitatively unchanged, but the 
depletion expands in space. This may reflect the 
observation from the movie 2 that the shape of 
the proton pulse is practically unchanged until 



tfl e w 2000. The proton pulse and the oxygen 
depletion are intertwined. 

The proton beam 1 and the electrons in our 
simulation have a significant mean speed in the 
reference frame of the proton phase space hole and 
they could provide the energy for the growth of the 
proton pulse. The movie 1 reveals that the proton 
phase space hole gives rise to vortices in the elec- 
tron phase space distribution and the electron flow 
is turbulent, which excerts a drag on the proton 
phase space hole. A coupling between the pro- 
ton phase space hole and the streaming electrons 
is established, which can also trigg er or accelerat e 
the growth of the phase space hole (|Dupredll986h . 
The simulation further evidences that the proton 
pulse grows at the expense of the kinetic energy of 
the proton beam 1. 

4. Discussion 

This work has examined secondary instabili- 
ties trigge r ed by the Buneman insta bility (BI) 
dBuneman 1958 ) with a PIC code ( Eastwood 
19911 ). One spatial x-direction has been resolved, 
which implies that the wave spectrum driven by 
the BI is monodirectional. This limits the realism 
of the simulation but it firstly allows us to resolve 
all relevant spatio-temporal scales and, secondly, 
it facilitates the identification of the relevant phys- 
ical processes. The large box has ensured that the 
periodic boundary conditions have not affected 
the results. We have considered two unmagne- 
tized, interpenetrating, uniform and unbounded 
beams that have been composed of three parti- 
cle species. The beam 1 has contained protons 
and electrons and the beam 2 has been formed 
by protons and by oxygen with a single negative 
charge. Each beam has been charge- and current 
neutral. The relative speed between both beams 
has been set to 0.4c and the initial temperature of 
all species has been selected to be 9.1 keV. 

The BI has saturated by t he trapping of elec- 
trons (jRoberts fc Berkl Il967l) . which groups the 
electrons into phase space holes that give rise to 
charge density modulations. The latter can be- 
have as quasi-particles and are thus referred to 
as plasmons. We could match the mean (group) 
speed of these plasmons with that of the protons 
of beam 2 by selecting negatively charged oxygen 
rather than electrons. Such a resonance facili- 
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tates the coupling of the p lasmon energy into the 
ion a coustic wave mode ( Mendonca fe Bingham! 
2002h . 

Initially the equilibrium between the elec- 
tric fields and the elect rons has been preserved 
( Korn fc Schamefll996allbl ) and the ion beams had 
lost only a small fraction of their energy; the elec- 
trons moved consequently only at mildly relativis- 
tic speeds. The thermal anisotropy of the electron 
distribution woul d trigger an ins tability similar to 
that proposed by Weibel ( 19591 ). resulting in the 
growth of magnetic fields and multi-dimensional 
field structures that cannot be represented by a 
ID PIC simulation. The Weibel instability may, 
however, not be important. The proton phase 
space holes form simultaneously with the collapse 
of the electron phase space holes that gives rise 
to the thermal anisotropy. The electromagnetic 
fields of the Weibel instability may not be strong 
enough to affect the interplay between the already 
present proton phase space holes and the oxygen 
and might be negligible. 

As the time progressed, two large density rar- 
efact ion pulses (I nfeld et al .ll989Hlnfeld fc Rowlands! 
1990) have developed in the protons of the beam 2. 
We have examined the larger one in more detail. 
The long time between the formation of the proton 
phase space hole and the pulse has suggested that 
it has been caused by the oxygen response. We 
have found that the proton density pulse grows, 
once the electric field of the proton phase space 
hole has displaced a substantial fraction of the 
ox ygen ions. A relat e d effe ct has been observed 
bv iEliasson fc Shuklal (|2004l ). There it was shown 
that an electron phase space hole was ejected by 
a density depletion in the background ions, which 
its electric fields excavated. Here, the consequence 
of the oxygen depiction is not the acceleration of 
the phase space hole but its growth. The electrons 
may also contrib ute to the gro wth of the proton 
phase space hole (jDupree Il986l) . 

The growth of the proton density pulse has re- 
sulted in an electron double layer that accelerated 
the electrons to ultrarelativistic speeds. To the 
best of our knowledge this is the first time that 
such an electron acceleration could be observed in 
PIC simulations of nonrelativistic beam instabili- 
ties. The electrons have reached a peak energy of 
50 MeV. This immense electron acceleration im- 
plies, that the ID PIC simulation geometry is not 



appropriate any more. Firstly, the electrons mov- 
ing at such speeds d rive the oblique mode insta- 
bility ( Bret A. 20061 ). rather than the two-stream 
instability. The oblique modes are excluded by 
the simulation geometry. Secondly, the extreme 
electron temperature anisotropy resulting from 
the thcrmalization of the ultrarelativistic electron 
beam would drive a very strong Weibel-type in- 
stability. The subsequent magnetic field growth 
would imply that electrons emit syn chrotron ra- 
diation on top of the bremsstrahlun g (jSchlickeiser 
l2003t iFleishman fc Toptvginll2007l). The bala nce 
of these processes (jSchlickeiser fc Lerche 2007 ) re- 
quires a three-dimensional geometry and a method 
to deal with the radiation, which is not possible 
with our PIC code. 

At the same time, this electron double layer 
constitutes a plasma structure that can easily ac- 
celerate the electrons to a speed, which allows 
them to enter the diffusive shock acceleration 
mechanism. This could be achieved with beam 
speeds that are not much high er than the flow 
speeds of supernova blast shells (jVink et al. 200G; 
iFransson et al.ll2002n . The oxygen ions we have 
modelled can also be found in such plasmas, al- 
though not in the high density an d in the ionisa- 



tion number we have assumed here ( Kat suda et al 



2008). Our negative charge is more representa- 



tive of dust, which we als o find in supernova blast 
shells (jMeikle et al.1 120061 ). Our findings can thus 
not be applied directly to a supernova scenario. 
However, the formation of the electron double 
layer in the presence of heavier ions clearly demon- 
strates the need to introduce these into PIC simu- 
lations, which so far only consider one ion speci es 
(jAmano fc Hoshinoll2007l ; IChapman et al.ll2005l) . 

Future simulation work must address the min- 
imum speed, at which the double layer can form. 
This has to be done in form of parametric sim- 
ulation studies, because it is not clear for which 
combinations of the beam speeds and tempera- 
tures th e double layers can form. Existing analytic 
studies (jVerheest fc Pillavi 120081 ) do not cover the 
case we consider here. The effects of the obliquely 
propagating modes driven by the initial beam in- 
stability have to be examined, as well as those 
due to the Weibel instability driven by the ther- 
mal anisotropy of the electron distribution after 
the initial BI has saturated. This requires at least 
two-dimensional PIC simulations. Finally, other 
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ion species and plasma compositions have to be 
examined to assess how easily double layers can 
form. This is important, bec ause the ion compo - 
sition is not easily obtained ( Ellison et al. Il99?t ) 
and not unique. 

The modelling of double layers and their conse- 
quences is also of interest for other astrophysical 
scenarios. Ion double layers, which are character- 
ized by ion distributions that resemble the electron 
distribution here, may regulate the ma ss inflow 



onto neutron stars ( Williams et al 



19861) and ac- 



celera te particles in accretion discs (M oslem et al 
20071 ). T hey may also contribute t o solar coronal 
ejections (jAlfven fc Carlquistill967[ ). 
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Fig. 1. — The linear growth rate f2s/f2 e is given 
as the height and as a function of the flow-aligned 
wavenumber Z x = k x Vb/VL e and one perpendic- 
ular Z z = k z Vb/£l e . The growth rate maximum 
ils/Qe = 0.049 is found at Z z = and Z x = 
0.59, which is close to the Jlsc/^e = 0.056 and 
k u Vb/£l e = 0.5 of the cold BI. 




Fig. 3. — The phase space distribution of the 
electrons (a) and of the protons of beam 2 (b) 
at tVl e = 145: The colour scale shows the 10- 
logarithm of the number of CPs. The electrons 
form phase space holes, which are coalescing. The 
horizontal line corresponds to the mean momen- 

1/2 

turn Vb/(1 — vl/c 2 ) of beam 2. The protons 
show velocity modulations that are spatially cor- 
related with the electron phase space holes. 
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Fig. 2. — The electrostatic energy density Ep and 
the kinetic energy densities Ei of the four particle 
species i: Panel a) shows Ep and the Ex of the 
electrons of beam 1, both normalized to E\(t = 0). 
Panel b) plots the E2 of the protons of beam 1, 
the E3 of the protons of beam 2 and the E4 of 
the oxygen of beam 2, which are all normalized to 
their respective initial value. 





Fig. 4. — The phase space distribution of the 
electrons (a) and of the protons of beam 2 (b) 
at t£l e = 10 3 : The colour scale shows the 10- 
logarithm of the number of CPs. The electrons 
display a plateau distribution for x/\ u > 10 and 
vortices for x/X u < 8. The plateau is split at 
x/X u 9. The proton beam shows tenuous fil- 
aments, which are protons trapped in the ion 
acoustic waves. A phase space hole is forming at 
x/X u ~ 9. 
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Fig. 5. — The normalized electrostatic field 
eE x /(Q e cm e ) (a), the proton density of beam 2 
(b) and the electron density (c) . The electric field 
shows a bipolar oscillation at x/X u ss 9, which 
modulates the densities of the protons and elec- 
trons. The time is tO e = 10 3 . 
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Fig. 7. — The electron distribution at tCl e = 2600. 
The phase space distribution in (a) demonstrates 
that the width of the double layer is still 5A n , but 
that the electrons reach now 50 MeV. The elec- 
tron phase space distribution is intergrated over 
the left interval x/X u < 50 in (a) and the number 
density is plotted in (b). Two beams are visible. 
The integral over x/X u > 66 in (c) shows a single 
electron beam. 





Fig. 6. — The phase space distributions at tfl e = 
2000 of the electrons (a), the protons of beam 2 
(b) and of beam 1 (c). The colour scale is the 
10-logarithm of the number of CPs. The electron 
distribution evidences a double layer in the inter- 
val 15 < x/X u < 20. The Double layer is main- 
tained by the proton pulse moving with the beam 
2, which is strong enough to modulate the protons 
of beam 1 . 
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Fig. 8. — The proton distribution at t£l e = 2600 
of beam 2 (a) and of beam 1 (b) . Panel (c) shows 
the normalized electric field eE x / (£l e cm e ) . A large 
electrostatic shock is forming at x/X u ~ 55, which 
involves the protons of both beams. An electro- 
static shock involving only the protons of beam 2 
is driven by the rapid expansion of the pulse at 
x/X u w 92. 
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Fig. 9. — The oxygen density at the times tO e = 
1000 (a), at m e = 1600 (b) and m e = 2000 (c) 
in the reference frame convected with the proton 
beam 2. The location of the original proton phase 
space hole is at (x — Vbt)/X u w 20. The density 
is integrated over 3 grid cells to reduce the noise. 
The density in (a) shows a density depletion of 
the order of the number density fluctuations. The 
density in (b) is modulated by about 50% and it 
has expanded in size in (c). The density varies 
approximately linearly over a limited interval in 
(b) and in (c). 
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